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We consider the flow of polarization current J = dP/dt produced by a homogeneous electric field 
£(t) or by rapidly varying some other parameter in the Hamiltonian of a solid. For an initially 
insulating system and a collisionless time evolution, the dynamic polarization P(t) is given by 
a nonadiabatic version of the King-Smith- Vanderbilt geometric-phase formula. This leads to a 
computationally convenient form for the Schrodinger equation where the electric field is described 
by a linear scalar potential handled on a discrete mesh in reciprocal space. Stationary solutions in 
sufficiently weak static fields are local minima of the energy functional of Nunes and Gonze. Such 
solutions only exist below a critical field that depends inversely on the density of k points. For 
higher fields they become long-lived resonances, which can be accessed dynamically by gradually 
increasing £ . As an illustration the dielectric function in the presence of a dc bias field is computed 
for a tight-binding model from the polarization response to a step- function discontinuity in £ (£), 
displaying the Franz-Keldysh effect. 



PACS numbers: 71.15.Qe, 78.20.Bh 
I. INTRODUCTION 

A very successful theoretical and computational frame- 
work was developed by King-Smith and Vanderbihpi for 
dealing within periodic boundary conditions with the 
macroscopic dielectric polarization of an insulator. The 
central result of the theory of bulk polarization (TBP) 
is an expression for the electronic contribution P which 
takes the form of a Berry's phased of the valence-band 
Bloch wave functions transported across the Brillouin 
zone (BZ). Alternatively, it can be recast in real space 
as the vector sum of the centers of charge of the valence- 
band Wannier functions. Practical prescriptions were de- 
vised for computing both the Berry's phasei and the 
Wannier functions, 3 which have become standard fea- 
tures of first-principles electronic structure codes. 

The measurable quantity accessed by the TBP is 
the change AP in macroscopic polarization induced by 
changing some parameter A in the electronic Hamilto- 
nian H(t) — H[X(t)]. The following assumptions were 
explicitly made in the original derivation 1 . (i) Adiabatic- 
ity: the change in X(t) is slow enough that the electrons 
remain in the instantaneous ground state of H(t), apart 
from small deviations proportional to dX/dt described by 
first-order adiabatic perturbation theory; (ii) the ground 
state of H(t) remains insulating at all times, separated 
from excited states by finite energy gaps; (iii) H(t) is 
lattice-periodic. The first two assumptions are related in 
that the size of the energy gap sets the scale for devia- 
tions from adiabaticity. 

A spatially homogeneous electric field necessarily vi- 
olates either (i) or (iii): if the field is introduced via 
a vector-potential A(t) = —c J £(t')dt', Hit) remains 
lattice-periodic but changes nonadiabatically, even for a 
static field; if instead a scalar potential term e£ (t) • f is 
used, H(t) is no longer lattice-periodic. Nevertheless, the 
TBP has been successfully applied to situations where 



electric fields are present ^^l&Sil^ but a rigorous justi- 
fication for doing so is still lacking. 

In this paper we reexamine the TBP and find that it 
can be generalized as follows. Assumption (i) can be 
dropped altogether. Assumption (ii) is only invoked at 
t = 0; the ensuing nonadiabatic dynamics may admix 
considerable amounts of excited states into the occupied 
subspace. Finally assumption (iii) can be relaxed to allow 
for a linear scalar potential to be present in addition to 
the periodic crystal potential. 

These generalizations extend the scope of the TBP 
to nonadiabatic polarization currents induced by time- 
dependent electric fields, or by other rapid changes in 
H{t) (e.g., the initial nonthermal ionic motion that ac- 
companies photoexcitation of the electrons by an in- 
tense laser pulse^i) The dynamical equations for the 
electrons that come out of this generalized TBP are 
derived and applied in the context of a tight-binding 
model. These equations are semiclassical (the electrons 
are treated quantum-mechanically, whereas the electric 
field is treated classically) and nonperturbative (electric 
fields of finite magnitude are allowed). 

We begin by considering in Section [D] some general 
properties of the coherent dynamics of Bloch electrons 
that are initially in an insulating state. They are used 
in Section ITTT1 to discuss the macroscopic current J(i), 
which is expressed as the rate of change of a dynamic 
polarization P(t) given by nonadiabatic versions of the 
King-Smith- Vanderbilt expressions. In Section IIVI we 
derive from this generalized TBP a numerically conve- 
nient form for the time-dependent Schrodinger equation 
(TDSE) in the scalar-potential gauge, discretized on a 
mesh of k points. Stable stationary solutions in static 
fields are discussed in Section They exist only below 
a critical field £ c which decreases with increasing fc-point 
density, and are local minima of the energy functional of 
Nunes and Gonze A prescription is given for com- 
puting them using an iterative diagonalization scheme. 
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In Section IVII we show numerically on a tight-binding 
model how the regime above the critical field can be ac- 
cessed dynamically, by gradually increasing the electric 
field beyond the critical value. We also compute the di- 
electric function of the same model in the presence of a 
static bias field, displaying the Franz-Keldysh effect. 



II. GENERAL PROPERTIES OF THE 
DYNAMICS 

A. Lattice-periodicity 

Here we expound in more detail an argument, sketched 
in Ref. la, that makes use of the one-particle density ma- 
trix to handle the presence of electric fields in a well- 
controlled fashionii^ We say that the one-particle density 
matrix n(r, r') = (r|n|r') is lattice-periodic if 

n(r,r') = n(r + R,r' + R), (1) 

where R is a lattice vector. In particular, this implies 
periodicity of the charge density. Suppose that is 
true at t = (e.g., the electrons are in the ground state 
of the crystal Hamiltonian H°(t = 0)). At that time a 
homogeneous electric field is turned on, which may sub- 
sequently have an arbitrarily strong and rapid variation; 
H°(t) may also undergo arbitrarily rapid variations (but 
must remain periodic). The full Hamiltonian in the scalar 
potential gauge is 

H(t)=H°(t)+H £ (t), (2) 

where H £ {t) = e£(i) ■ r describes the electric field in the 
dipole approximation, and — e is the electron charge. 

Let us show that in the absence of scattering the 
lattice-periodicity of n(r, r') is preserved at all later 
times. It suffices to establish that n(r, r') = n(r + 
R, r' + R) . The density matrix evolves according to 
ihdn/dt = [H,n], or, in the position representation, 

ihh(r,r') = J [H(r, x)n(x, r') - n(r, x)i?"(x, r')]dx. (3) 

(When left unspecified, the domain of integration over 
spatial coordinates is understood to be the entire space.) 
For clarity we consider the effect of H° and H sepa- 
rately. The H° term yields 

z^n(r + R,r' + R) = J [H°(r + R, x)n(x, r' + R) 

- n(r + R,x) J ff°(x,r' + R)]dx(.4) 

Making the change of variables x' = x — R and invoking 
the lattice-periodicity of H° and h, we find n(r + R, r' + 
R) = n(r, r'). Using r(r,r') = (r|f|r') = rS(r — r') the 



contribution from H £ is seen to have the same property: 

iS,7i(r + R,r' + R) = e£ ■ (r + R) n(r + R, r' + R) 
- e£ ■ (r' + R) n(r + R, r' + R) 
= e£ ■ (r — r') n(r, r') = ifi n(r, r'). 

(5) 

Hence n(r,r') remains lattice-periodic under the action 
of the full Hamiltonian @. This was to be expected, 
since in the vector potential gauge the Hamiltonian is 
periodic^ The purpose of this exercise was to show 
explicitly how this result comes about in the scalar- 
potential gauge, where the nonperiodicity of H has been 
a source of some confusion regarding this issue. 

B. Wannier-representability 

The previous result on the conservation of lattice- 
periodicity is valid for both metals and insulators. In 
what follows we shall specialize to the case where the 
system is initially in an insulating state, in which case a 
stronger statement can be made regarding the nature of 
the states at t > 0. 

We will assume the absence of spin degeneracy 
throughout, so that states are singly-occupied. In terms 
of the valence Bloch eigenstates of H°(t = 0), the initial 
density matrix is 

M „ 

n(r,r';t = 0) = n B 1 Y, / dkVkn(r)Vk„(r'), (6) 

n=l 

where the integral is over the BZ of volume Qb = 
(2tt) 3 /v, and M is the number of filled bands. (Clearly, 
such a density matrix is lattice-periodic. Its idempotency 
can be checked using Eq. HAIL ) We shall prove that, as 
the density matrix evolves in time according to 

ihh{r,r';t) = (r\[H° + H £ , n]\r'), (7) 

it can still be expressed in the same form, 

M 

n{v,v'-t) = n- 1 Y, / dk0 kn (r,i)0* n (r',i). (8) 

n=l •* 

Although at t > the occupied states 0kn(r,£) may 
depart significantly from the valence states of H a (t), 
they remain orthonormal and Bloch-like: (j)k n (r,t) = 
e lU r v Wn (r,t), with w kn (r + R, t) = v kn (r,t). 

The cell-periodic states Vk n (r,t) are the central ob- 
jects in our formalism. For discussion purposes only, 
let us expand them in the set of eigenstates Mkm(r,i) — 

e _lk ' r, 0km(r, t) of the cell-periodic Hamiltonian H^(t) = 
e -ik.?£o^ykf. 

OO 

\Vkn(t)) = ^ C ^,nm{t) \u km (t)}. (9) 
m— 1 
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Individual eigenstates will in general have fractional oc- 
cupations < n km = J2n=l \ c Knrn\ 2 < 1 at t > 0, but 

the total population = J2m=i "km is the same for ev- 
ery k and equals the number of filled bands at t = 0. 
This is intuitively clear, since a spatially homogeneous 
electric field causes vertical transitions in fc-space which 
amount to a redistribution of the electron population 
among states with equal k; the same is true for the tran- 
sitions induced by varying the lattice-periodic H°(t). 

We will justify Eq. (JSJ by deriving a dynamics for the 
\ v kn) that insures that Eq. (JSJ provides a solution to (7J|. 
Since there is a gauge freedom 



M 



(10) 



([/k is a k-dependent unitary M x M matrix) in the def- 
inition of the |i>kn}^ the evolution equation for them is 
not unique. We require only that the |t>kn) should yield 
the correct dynamics for the gauge-invariant density ma- 
trix, Eq. 0, and we will look for the simplest solution 
that achieves this goal. 

By hypothesis, at time t n(r,r') takes the form 



n(r, r') 



M 



[Wkn(r)w kn (r') +Vk„(r)< n (r')]. (11) 



As in the previous subsection, we consider the contri- 
butions from H° and H in Eq. Q separately. The 
former is captured by iTi\v-^ n ) — H^\v-^ n ). To deal with 

H we resort to manipulations familiar from the crystal- 
momentum representation (CMR)i^ (but with the crucial 
difference that in the CMR those manipulations are ap- 
plied to the \uk n ), not to the |ukn))- We first observe 
that 

(r| [f, n] |r') = (r — r') n(r, r') 

M 



will satisfy Eq. J7J). The time-independent version was 
introduced as an ansatz in Ref. The equivalence of 
Eq. H14JI to other forms in the literature is established in 
Appendix ^ 

If at time t the M states \vk n ) at every k are lattice- 
periodic and orthonormal, the dynamics dictated by 
Eq. itH|l preserves those properties, i.e., Ukn(r + R) = 
Wkn(r) and d(vk n \Vk m ) / dt — 0- This can be seen as fol- 
lows. Starting from 

ihv kn (r+R) = J Fg(r+R,x)«kn(x)tix+ie£-9k«kn(r+R), 

(15) 

making the change of variables x' = x — R, and invoking 
the assumed lattice-periodicity of both H° and Ukn(r), 
the right-hand-side becomes ihvk n (r). As for orthonor- 
mality, Eq. ifHj) yields^ 

d e 

^(Wknbkm) = ■ dk{Vkn\Vkm) ■ (16) 



Since by hypothesis 



{Vkn \V kr, 



«kr 1 ( r )«km(r)c?r 



(17) 



where the integral is over a unit cell, the right-hand-side 
of Eq. I)lri|l vanishes. This completes the proof of Eq. JBJ). 

Two assumptions were made in the above derivation. 
The first is that the states |«kn) vary smoothly with k, 
so that /c-derivatives exist; we will come back to this 
point in Sec. IIII Bl The second is that the dynamics 
is scattering-free. Note that Eq. (|16[) is closely related to 
the collisionless Boltzmann equation; incoherent scatter- 
ing would destroy the constancy of the total population 
nk by inducing transitions between different k points. 

Having established that the occupied manifold is 
spanned by M Bloch-like states at each k, we now 
transform them into Wannier-like states (r|WR„(f)} = 
W n (r — R, t) in the usual way, 



" l r m 

n—1 I J 

rn—1 



dke- ikR C/ k .. 



»(*) \fam(t)), 



Integrating by parts and noting that in a periodic gauge 
(0k+G,7! = </^kn) the boundary term vanishes, we obtain 

M . 

(r\[H £ ,h]\r') =n B 1 Y, / dke* {r - r,) ie£- 

n=l •* 

"(9k U k n (r))< n (r')+^k„(r)(5k Ukn (r'))](13) 

Comparing with Eqs. J7J) and (|11() we arrive at ih\v^ n ) = 
ie£ ■ <9k|«kn}- The effect of H £ thus takes the form of a 
fc-derivative, and the combined effect of H° and H £ is 



where a periodic gauge is assumed and we have inserted 
a unitary rotation l|10|) among the occupied states. The 
assumption that by a judicious choice of the matrices 
U-k{t) the Bloch-like states can be made to vary smoothly 
with k is equivalent to the assumption that the Wannier- 
like states can be chosen to be well localized; 3 
The density matrix © can now be recast as 



M 



i(r, v'-t) = £ W *^ r > *) W ^ r '> *)■ 



(19) 



n=l R 



ih\vkn) = (Bk + • dk)\v kn ). 



(14) 



This is our version of the TDSE for Bloch electrons in the 
scalar potential gauge, constructed in order that Eq. © 



We will term Wannier-representable (WR) a state whose 
density matrix is of this form. An insulating ground state 
is WR, while a metallic state is not. We have established 
that under the Hamiltonian @ and in the absence of 
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scattering, an initially insulating system remains WR, or 
"insulating-like" , even if at some later time the ground 
state of H° (t) becomes metallic^ Unlike a true insulat- 
ing ground state, or a stationary field-polarized stated 
a dynamic WR state will in general break time-reversal 
symmetry and carry a macroscopic current. This is the 
subject of the next Section. 



III. DYNAMIC POLARIZATION AND 
CURRENT 

A. Derivation 

Our aim in this Section is to show that a WR state 
carries a current that can be expressed as the rate of 
change of a polarization per unit volume, 



J(t) = 



dP(t) 
dt 



(20) 



where P(£) is given in a periodic gauge by 



P a (t) 



(2tt) 



/ dk{v kn (t)\d ka v kn (t)) (21) 



(a is a cartesian direction) or, equivalently, by 



M 



P(t)=-- J2 / r\W n (r,t)\ 2 dv 



(22) 



Eqs. (j2*T|) - (|2*2*j) are identical to the King-Smith- Vanderbilt 
expressions appropriate for the adiabatic regime and £ = 
Ofi except that in 121JI the valence-band eigenstates |itkn) 
have been replaced by the instantaneous solutions of the 
TDSE JEJ, and the W n (r, t) in Eq. JUJ) are the Wannier 
states corresponding to the v kn (r,t). Eq. (|2~TJl can be 
interpreted as a nonadiabatic geometric phased 

As in the adiabatic case, P(i) is invariant under the 
transformation only up to a "quantum of polariza- 
tion" (e/u)R. Naturally, this gauge indeterminacy does 
not affect the measurable J(t). The total change in bulk 
polarization in a time interval [0, T] is also well-defined 

as the integrated current: AP = L 3(t) dt. It can be de- 
termined, apart from an integer multiple of the quantum, 
by evaluating P(t) at the endpoints: AP = P(T)-P(0). 
In practice the remaining indeterminacy can be removed 
in the manner described in Ref. 01 by evaluating P(t) 
with sufficient frequency during that interval. 

To establish Eqs. we first evaluate dP/dt by 

taking the time derivative of Eq. I|21|l and obtain, after 
an integration by parts, 



dP a 
dt 



M 



^2 dk[(v kn \d ka v kn ) - ex.]. (23) 



Inserting the TDSE, Eq. I|14f> . we note that the contribu- 
tion arising from the second term therein, which explic- 
itly involves £, may be written as 



M . 

£X>y dk O k ) ( 24 ) 



^ (27r) a ft ^ 
where 

^L^ ( k ) = i[(dk a Vkn\dkpVkn) ~ (dk0Vkn\dk a Vkn)] • (25) 

This takes the form of a (nonadiabatic) Berry 
curvature^ Using Stokes' theorem, its volume integral 
can be turned into a surface integral around the edges of 
the BZ of the Berry connection Ak ««, where 



A k.,nn = i(vkm|<9fe Q V k n)- 



(26) 



Such an integral vanishes in a periodic gauge, so that 
J a = 0. The remaining contribution, arising from the 
insertion of the first term of Eq. I|14|l into Eq. 1)23(1. then 
gives 

dP M r 

^f = -^E J dk[(v kn \HZ\d kc ,v kn )+c.c.]. (27) 



(2*)'.. „=i 
On the other hand, the current is 



J a = --Tr e (nv a ). 
v 

Here Tr c denotes the trace per unit cell, 
Tr c (0) = ^ f 0(r,r)dr, 



(28) 



(29) 



where N is the (formally infinite) number of real-space 
cells in the system. The velocity operator is defined as 



1 r 

v a = — [r a ,H\. 
in 



(30) 



Inserting the Hamiltonian J5J and using [f a ,H £ ] = 0, 



v a = -[f a ,H ]. 
in 



(31) 



In the position representation we find, combining JSj, 
(|28|l and l|31|l . invoking the lattice-periodicity of the in- 
tegrand to replace (1/iV) J dr by J dr, and inserting the 
identity i = Jdr'|r')(r'|, 



J\ a 



~f^^ I dk I v dr I dr ' Uk «( r ') Wk «( r ) 



(2 



(32) 



Integrating by parts in k a (the boundary term vanishes 
in a periodic gauge), and using 



tf£(r',r) = e- ,;k -( r '^ff°(r',r), 



(33) 
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J a reduces to exactly the same expression appearing on 
the right-hand side of Eq. H27|) . This completes the proof 
of Eqs. f° r states evolving under the Hamil- 

tonian (J2J. 

We note in passing that the integral on the right-hand 
side of Eq. (|27J) can be recast as 



dk[d ka (vwn\H°\v kn ) - (vwn\(d ka H£)\vwn)]- (34) 



The first term vanishes in a periodic gauge, leading to 
the more familiar-looking form 



M „ 

/ dk(v kn \v a (k)\v kn ), (35) 

( 27r ) n=l J 



where v a (k) = {l/h)(dk a f$)M 

The above derivations (and indeed all the results in this 
paper) remain valid for nonlocal pseudopotentials such as 
those used in ab initio calculations, since the definition 
of the velocity as the commutator l|30|l remains valid for 
such pseudopotentials. 



B. Discussion 

It is remarkable that a knowledge of the wave func- 
tions at t — and t = T is sufficient to infer, to within a 
factor of (e/v)H, the net amount of current that flowed 
through the bulk in the intervening time. This is a di- 
rect consequence of representability by localized Wannicr 
functions, that is, of the insulating-like character of the 
many-electron system. For such systems the integral in 
Eq. H22fl can be evaluated, and it becomes possible to 
track the time evolution of the electronic center of mass, 
i.e., of P. Indeed, the center of mass can be meaning- 
fully defined within periodic boundary conditions only 
for many-electron states that are localized in the man- 
ner of insulating statesi^fliSi Under these conditions, the 
history of the coherent current flow is contained (mod- 
ulo the quantum of polarization) in the initial and final 
wave functions, related by the time evolution operator 

exp[-(i/h) £ H(t)dt]. 

This result was previously established for adiabatic 
charge flow, » under the assumption that the ground state 
is separated from excited states by finite energy gaps ev- 
erywhere in the BZ. In nonadiabatic situations the oc- 
cupied manifold acquires a significant excited-state ad- 
mixture, so that it becomes impossible to identify an en- 
ergy gap. Instead, underlying the derivation in Sec. lIII A1 
is a weaker assumption, namely, that the many-electron 
state has a localized nature, as reflected by the ability 
to construct, via Eq. 118(1 . Wannier functions having a 
finite localization lengthiSLS^ (Numerical calculations of 
the localization length will be presented in Sec. IVII ) For 
instance, when taking fc-derivatives, we assumed a "dif- 
fcrcntiable gauge" for the \vk n ). This is only possible if 
the character of the electronic manifold changes slowly 



with k, which is precisely what is measured by the lo- 
calization length. 3 - These observations are in line with 
Kohn's viewpoint that the defining feature of the insulat- 
ing state is wave function localization, not the existence 
of an energy gapi2£ 



IV. DYNAMICAL EQUATIONS 

Having found the Berry-phase formula l|21|) for the dy- 
namic polarization in the presence of a field £(t), let us 
now use it to obtain computationally tractable dynam- 
ical equations under the Hamiltonian (J2J. The starting 
point is the observation that the dipole term H £ con- 
tributes — vP(t) ■ £{t) to the energy per unit cell. An 
energy functional valid for periodic boundary conditions 
is then obtained by expressing P(t) via the TBP formu- 
las. This program was previously carried out for insu- 
lators in static fields where stationary states were 
computed by minimizing that energy functional after ap- 
plying a regularization procedure (truncation of the Wan- 
nier functions in real space or discretization of fc-space). 
Our strategy for the time-dependent problem is to im- 
pose stationarity on the corresponding action functional. 
Following Refs.Qla we adopt here a fc-space formulation, 
which is particularly well-suited for numerical work. Spe- 
cial emphasis will be put on the discrete-fc case since this 
is the relevant one for numerical implementations. 



A. Continuum-fc case 

In the continuum-fc limit the TDSE may be formally 
obtained from a Lagrangian density C(k) such that the 
Lagrangian per unit cell is L = fl^ 1 J dkC(k). For WR 
states under the Hamiltonian (0 we have 



M 



£(k) =ih^2 (v kn \vkn) - E(k), 



(36) 



71=1 



where 



A I 



E(k) = J2(vwn\H£ + ie£-du\v kn ). (37) 



n=l 



Using Eq. (|21|l and defining the zero-field energy func- 
tional 

M 

E°=n^Y^ dk(v* n \H°\v kn ), (38) 

n=l 

one finds the total energy functional 

E = Qg 1 J dkE(k) = E° - vP ■£. (39) 

The Euler-Lagrange equation 23 

d SC d 5C SC 



dt (flukwl dkiSdkVknl (Svu, r , 



(40) 
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then leads to the dynamical equation (|14|l . 

As already mentioned, the choice of dynamical equa- 
tion for the \vk n ) is not unique. An alternative to Eq. 114(1 
is 



iH\vkn) = (#k + ie £ ■ <9k)|vkri}- 

The bare derivative dk has been replaced by 

A I 



^ ^ Ak,mn |^km) (^k7i 



(41) 



(42) 



where Ak, m n is given by Eq. (J2EJ). The operator <9k is 
a multiband version of the covariant derivative^ and is 
discussed further in Appendix [BJ Although the field- 
coupling term in Eq. (|41|l is no longer a scalar potential 
term in the strict sense, we will continue to view it as 
such in a generalized sense. 

Eq. (|41(l preserves the orthonormality of the |fkn) and 
generates the correct dynamics for the density matrix. 
(These properties rely on (^4£)^ = ^k' which follows from 
dk a (vkn\vkm) = 0^) The latter is most easily seen from 
the dynamics of the projector 



A I 



Pk = ^2 \vkn)(Vkn\, 



(43) 



n=l 



which completely specifies the occupied subspace at k 
while being insensitive to unitary rotations inside that 
subspace. After some algebra, it can be shown that while 
the individual \vk n ) behave differently under Eqs. itT4*|) 
and l|4*T)> . Pk stays the same. 

An advantage of introducing Eq. i|41[) in place of 
Eq. HU) is that, upon the discretization of /c-space, the 
former leads to an evolution equation at point k that is 
gauge-covariant (in the sense of transforming in the obvi- 
ous way under unitary rotations among occupied states 
at k and being invariant under such rotations at neigh- 
boring points k'), as will become clear in the next section. 



B. Discrete-fc case 

This is the relevant case for numerical work. The La- 
jrangian for a uniform mesh of N points in the BZ is 



it) 



A I 



N 

n=i k 

where E is the energy in an electric field, 
E = E° - vS ■ P, 

with 

M 

N 



1 M 



(44) 



(45) 



(46) 



and a discretized expression for P to be given shortly. 
Applying the Lagrangian equations of motion 23 



d 5L 



SL 



dt (5v kn \ (dvkn\ 







yields 



Writing 



in j t ^ n) = - Nv£ ■ wL 



£l>(p-b.) 

i=l 



(47) 



(48) 



(49) 



where and are the direct and reciprocal lattice vec- 
tors respectively, and defining 



v P • b, = -er\ ; 
the last term in Eq. 148|) becomes 



(50) 



(51) 



According to the Berry-phase theory of polarization^ 
Ti is the string-averaged discretized geometric phase 
along the direction, 



1 

"N? 



E 

i=i 



N;' 



Imln Yl det5(k^ ) ,k^ 1 ). (52) 

j=0 



Here S mn (k, k') = (^kml^k'n) is the M x M overlap ma- 
trix, is the number of strings along bi , each contain- 
ing Ar| points k^ = k^ 



jAki, k^ is a point on the 



(b 2 , b 3 ) plane labeled by I, and Aki = bi/Aj 1 . Eqs. 
and (|52|l provide the discretized version of the nonadia- 
batic Berry-phase polarization, Eq. (|21|l . (A discrete-fc 
formula for the macroscopic current 3(t) is given in Ap- 
pendix [DJ) As in the continuum case, a periodic gauge 
is assumed. 

A compact expression for STi / (Svkn | is derived in Ap- 
pendix EH using the following notation. Let kicr = 
k + a Ak.i, where a = ±1. The overlap matrix becomes 
,mn — ( v km\vkia,n) ■ Next we define 



At 



Inkier, n) — ^ ' (^kicr 



(53) 



which is a "dual" of \vk n ) in the space of the \v)'s at the 
neighboring point kicr, since 



(54) 



n=l k 



The \vkia,n) are gauge-covariant in the sense that (i) they 
are invariant under unitary rotations among the \vkia,n) 
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at any neighboring point Via, and (ii) they transform un- 
der unitary rotations among the |t>kn) m the same man- 
ner as the |«kn) themselves, i.e., 

M 

|«kiff,n) — * 53 Uk,mn\vkia,m)- (55) 
m— 1 

Then it is shown in Appendix IC 21 that 

Wt\ = WEr a ^ itr ' n) ' (56) 
Combining Eqs. (|5T)l . and (JSHJl, and defining 

3 

h«k«) = j- 53 iV] '(5 • a*) 53 crlwkicr,™}, (57) 

1 = 1 (7=±1 

the dynamical equation becomes 

ift^_\vk n ) = Hk\v kn ) + \wkn}- (58) 

Eq. (|58|l is a discretized version of Eq. (|41|) . i.e., 

|w k „) ~ ie£ • 9 k |wkn}- (59) 

This is connected with the fact that the duals provide 
a natural framework for writing a finite-difference repre- 
sentation of d^lvkn}^ In our notation, 

). (60) 

1 <7=±1 

Both dynamical equations (|14|> and l|41|) lead to 
rf^knl^km}/^ = f° r WR manifolds, so that the time 
evolution of the individual states \vkn) is unitary. This 
property is preserved in the discretized form (|58|l . since 
(fkn|wkm) = 0. In order to take advantage of certain 
unitary integration algorithms, it is useful to recast the 
term \wk n ) on the right-hand-side as an hermitian op- 
erator acting on |^kn)- For that purpose let us define 



so that QkWk = Wk (where Qk = 1 — Pk) and therefore 
w^lvkri) = 0. We have thus achieved our goal: Eq. 
now takes the canonical form of a TDSE, 

■ih-^lvkn) = f k \v kn ), (63) 

with an hermitian operator on the right-hand side: 

fk(£) = H^ + w k (£)+wt(e). (64) 

We remarked previously that in the continuum-fc limit 
Eq. I|58(l reduces to Eq. 141|l . The corresponding analysis 
for Eqs. ftE ty -Ip p is left to Appendix El 

The operator Wk appearing in Eq. I|64l) is defined via 
Eqs. (J53J, (EU, and (J^J). It should be emphasized that 
it depends explicitly on the occupied states at k and Via. 
In particular, even when and £ are time- independent, 
if the occupied manifolds at k and Via are changing over 
time, so is T k - However, Tk remains invariant under uni- 
tary rotations at fc-points k and Via. Hence the resulting 
dynamics of the occupied manifold (Eq. (|66|l below) has 
the essential property of being insensitive to the gauge 
arbitrariness that is always present in numerical simula- 
tions. 

Note that when the Lagrangian procedure was applied 
in Sec. IIV Al to the continuum-fc problem, we arrived at 
Eq. which contains dk\vkn)- When the same was 

done after discretization, the resulting dynamical equa- 
tion contained instead <9k| u kn}- The reason is that the 
gradient of the discretized Berry's phase, Eq. is 
by construction orthogonal to the occupied subspace at 
k, whereas the corresponding continuum-fc term used in 
Sec. lIV Al was not. Had we orthogonalized that term, we 
would have obtained Qkdk\vkn) instead of dk\vk n ), which 
is equivalent to <9k|^kn) (see Appendix FBI . 



C. Numerical time integration 

In the applications of Sec. I VII we use the algorithm^ 



Vkicr,n){vkn\ 



(61) 



which converts an occupied state at k into its dual at Via 
and is invariant under gauge transformations (i.e., under 
independent unitary rotations among occupied states at 
both k and kicr). It follows that the operator 



Wk(£) 



47T 



53 Nf(£ -aOj^crfW 



(62) 



turns |i>kn) mto \wk n ), which is the property we seek. 
Lastly, for the purpose of acting on \vkn) the nonhermi- 
tian Wk can be replaced by Wk + w^ since PkPkia = Pk, 



\vkn(t + At)) = 



1 - ih(At/2)fk(t) 
l + ih(At/2)fk(t) 



Vkn(t)) (65) 



to perform the time integration. Note that in order to 
use this algorithm it was necessary to invoke the form 
(EH of the TDSE. The hermiticity of f k guarantees that 
the time evolution is strictly unitary for any value of 
At. Since the system under study in Sec. IVII is a tight- 
binding model with only three basis orbitals, the matrix 
inversion is very inexpensive. The same algorithm has 
been successfully used to perform self-consistent time- 
dependent density-functional calculations of the optical 
properties of atomic clusters using localized orbitals as 
a basis set4^ For calculations with large basis sets (e.g., 
plane-waves) more efficient algorithms are available 
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Owing to the hermiticity of Tk, the projector 143|) 
obeys 

Hence, a variation of the above approach would be to 
replace Tk by 

% = Qkf k + T k Q k , (67) 

which is also hermitian. Because [7k, = [Tk, this 
choice does not change the dynamics of the occupied sub- 
space Pk, but it does change the dynamics of the individ- 
ual states |i>kn)- In fact, 7k generates a parallel transport 
evolution characterized by (vkn|i>kn) = 0, thus discard- 
ing the "irrelevant" part of the dynamics associated with 
phase factors and unitary rotations inside the occupied 
subspace. We have found empirically, however, that the 
use of 7k in Eq. I|tj5[) appears to result in a less stable 
numerical time evolution, and we have therefore chosen 
to retain the original Tk dynamics in our practical imple- 
mentation. 

D. Discussion 

It may seem surprising that a linear potential can be 
accommodated in a theoretical description of a periodic 
bulk system. A commonly-held viewpoint is that a linear 
potential can be implemented within periodic boundary 
conditions only for the case of a finite system (molecule or 
cluster) in a supercell, in which case it becomes possible 
to introduce a sawtooth potential as long as its disconti- 
nuity is located in a region of negligible electron density. 
To the contrary, Eqs. (fT4*|) . (|4*T)) . and (|55|> demonstrate 
that it is perfectly permissible to insist on the usual pe- 
riodic boundary conditions on the wave functions while 
allowing for nonperiodicity of the potential. This can be 
done because the potential takes the special form of a 
sum of spatially periodic and linear contributions, rel- 
evant to a crystal in a homogeneous electric field. As 
shown in Sec. Ill Al the action of the nonperiodic Hamil- 
tonian, Eq. @, then preserves the lattice periodicity of 
the density matrix, which can therefore be represented 
by periodic wave functions. 

Incidentally, we note that a sawtooth operator of sorts 
is "hiding" behind the TBP formulas. The Berry-phase 
polarization has been recast as the expectation value of a 
properly-defined center-of-mass position operator of the 
many-electron periodic system^i That operator, intro- 
duced by Kohn, 20 is a sawtooth, not in real space but 
in the configuration space of the many-body wave func- 
tion. It can only be constructed for wave functions hav- 
ing a certain disconnectedness (localization) property in 
configuration space characteristic of the insulating state. 
This observation is closely related to the discussion in 

Sec. lrnBl 

Finally, we mention that an alternative approach 
for introducing a linear potential into a periodic solid 



is via the crystal momentum representation (CMR) 
formalism^ This approach is summarized in Ap- 
pendix A, where the connection with our formalism is 
established. The CMR dynamical equations appear to 
be less convenient for computational work. However, the 
advantages of the present formulation came at the ex- 
pense of generality, since our equations are restricted to 
the scattering-free dynamics of initially insulating sys- 
tems. 



V. STABLE STATIONARY SOLUTIONS 

A. Formulation 

Let us try to find, for a constant £ 0, solutions 
to Eq. i|t)3|) for which the occupied manifold remains 
unchanged over time. A natural guess is the manifold 
spanned by M eigenstates of Tk at each k, 

Tkl^kn) = -Ekn(£)Kn}- (68) 

Since Tk depends on the occupied states at the neigh- 
boring /c-points, Eq. I|68|l must be solved self-consistently 
among all k. If a solution exists, the corresponding Tk 
and i\ commute and, according to Eq. (|S5|> . dP^/dt = 0, 
i.e., the solution is stationary. 

We are now ready to make contact with Refs. I Slit 
where the energy functional E of Nunes and Gonze, 7 
Eq. (|45|) . was minimized at fixed £. A stationary point 
of E has zero gradient: \G kn ) = SE/(Svkn\ — 0, where 
the functional derivative is taken in such a way that the 
gradient is orthogonal to the occupied space. In Ap- 
pendix 1^^] it is shown that 

|Gk„> = {l/N)Q k f k \vk n }, (69) 

so that solutions of Eq. (|68|) obey |Gk n ) = 0. Thus, 
stationary solutions of the dynamical equation are sta- 
tionary points of E. 

A Hessian stability analysis 30 shows that a necessary 
condition for a stationary point of £ to be a minimum 
is that the M lowest-lying eigenstates of Tk are chosen. 
Since doing so at £ = yields the ground state, at fi- 
nite £ that procedure yields a state that is is adiabat- 
ically connected to it by slowly ramping up the field, 
keeping the system in a minimum of E. Such "polarized 
manifolds" have been discussed previously in a pertur- 
bative framework^SLSiiiS treating k as a continuous vari- 
able. In that limit the electric field perturbation becomes 
singular. That is, even an arbitrary small field induces 
a current via Zener tunneling to higher bands, and the 
polarized manifolds are not stationary, but rather, are 
long-lived resonances. In other words, an infinite crystal 
in the presence of a static electric field does not have a 
ground state. This is reflected in E loosing its minima 
as soon as £ departs from zero. 

Instead, for a discrete mesh of fc-points, arguments can 
be giveni& suggesting that E loses its minima only when 
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£ exceeds a critical value £ C (N) that decreases as the 
number TV of k points increases; this is supported by nu- 
merical calculations. 8-9 It follows from the preceding dis- 
cussion that the minima of E below £ C (N) are stable sta- 
tionary solutions of the dynamical equation. Conversely, 
above £ C (N) there are no such solutions. 

These two regimes - below and above the critical 
field - will be explored numerically in Sec. I VI CI via 
time-dependent calculations. If one stays below £ C (N), 
the stationary solutions can be computed using time- 
independent methods, such as the diagonalization algo- 
rithm described next or the minimization methods of 
Refs.Hl 



B. Diagonalization algorithm 

We have in Eq. (|68(l the basis for an algebraic method 
of computing stationary states at finite £ on a uniform 
fc-point mesh, for \£\ < £ C (N): loop over the k points; 
for each one select the M eigenstates of with the low- 
est eigenvalues; iterate until the procedure converges at 
all k and the occupied subspace stabilizes (this will only 
happen below £ c ). Even in a tight- binding model with- 
out charge self-consistency, the set of equations H68JI has 
to be solved self-consistently throughout the BZ, since 
the operators Tk couple neighboring k points via their 
dependence on the |ukn}- One may choose to update Tk 
either inside or outside the loop over k; the latter option 
renders the algorithm parallclizable over k points. 

We have tested this scheme on the tight-binding model 
of Sec. I VII and confirm that it produces the same state 
as a direct steepest-descent or conjugate-gradients mini- 
mization of the functional E& This algorithm may be es- 
pecially suited for implementation in certain total-energy 
codes that are based on iteratively diagonalizing the 
Kohn-Sham Hamiltonian expanded in a small basis set 
of local orbitals^ 



then find 

M . 

n^ 1 / dkE kn (£) = E°{£) >E°(£ = 0), (71) 

n=l 

where E°(£) is the zero-field energy functional (|38|) eval- 
uated at the field-polarized stationary state, and the 
inequality follows from the variational principle. The 
same properties hold for the eigenvalues of the discretized 
form l|68[l . which can be obtained by diagonalizing 
inside the occupied manifold. We have here the inter- 
esting situation that a minimum of the total energy E 
can be obtained by solving the eigenvalue equations l|68() 
whose eigenvalues, summed over n and k give instead 
the zero-field contribution E°. This can be traced back 
to Eq. (|B2(1 , which expresses the "parallel-transport- like" 
nature of the covariant derivative. 

The above is to be compared with the time- 
independent version of Eq. (|14|l , 

+ ie£ ■ d k )K„) = £kn(£)K»>- (72) 

Under a diagonal transformation \v' kn ) — > e l9km \v' kn ) its 
eigenvalues change as E^ n (£) — > E kn (£) — e£ ■ 9k#kn- 
The analog of Eq. l(7T|) is 

M 

fl^ 1 £ / dkEUe) = E\£) - v£ ■ P(£). (73) 

n=l 

The quantity on the right-hand side is now the total en- 
ergy E, which is invariant only modulo e£ ■ R. 

Although the individual eigenstates of Eq. I|7l)|l are 
in general different from those of Eq. (|72[1 . the self- 
consistent solutions for all k and n span the same space 
in both cases, i.e., they differ only by a gauge transforma- 
tion. It is then a matter of convenience to choose which 
of the two equations to solve in practice. Our particular 
approach is to discretize Eq. lf7U|) in a gauge-covariant 
manner, and then solve the resulting Eq. I|68|) . 



C. Discussion 

Eq. (|68|l is a discretization of the time-independent ver- 
sion of Eq. lj4"T|) : 

(Hi + ie£ • d k )Kn> = E kn {£)\v kn ). (70) 

An analysis of the eigenvalues of this equation will serve 
as a guide for discussing those of Eq. I|68|) . (For the 
present purposes we will assume that the continuum 
form H7(J[) has solutions for f ^ 0.) As a result of the 
properties of the covariant derivative (Appendix [BJi the 
Ek n (£) are invariant under diagonal gauge transforma- 
tions C/k,mn = e l9km ^m,n- Upon multiplying on the left 
by (vk n | the second term on the left-hand-side of Eq. I|70|) 
vanishes. Integrating over k and summing over n, we 



VI. NUMERICAL RESULTS 

A. Tight-binding model 

We have applied our scheme to the one-dimensional 
tight-binding model of Ref. 0, a three-band Hamiltonian 
with three atoms per unit cell of length a = 1 and one 
orbital per atom, 

H°(a) - teWfe + h.c] }, (74) 

3 

with the site energy given by £3 m +;(a) = Acos(a — j3{). 
Here m is the cell index, I = { — 1, 0, 1} is the site index, 
and f3i = 2ttZ/3. Before the Berry-phase polarization 
can be computed 34 (or an electric field applied to the 
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FIG. 1: Energy dispersion of the tight-binding model for the 
choice of parameters t = 1, A = — 1, and a — 0. 

systemic), the position operator must be specified. Al- 
though this may be done without introducing additional 
parameters^ we adopt the simple prescription of Ref. 4: 
x = Xj&jCj, with Xj — j/3. In the results reported 
below we have set e = h = 1, t = 1, and A = — 1, and 
only the lowest band is filled (with single occupancy). 
Fig. ^ shows the band structure at zero field for a = 0. 



B. Sliding charge-density wave 

The Hamiltonian of Eq. H74|) is a simple model of a 
commensurate charge-density wave which slides by one 
period as the parameter a evolves adiabatically through 
2tt. It is easiest to see this by noting that in the space 
of parameters A x = A cos a and A y — A sin a, cycling 
a by 2ir corresponds to tracing a circle about the origin 
in the A^-Aj, plane. The system is insulating (i.e., a 
gap remains open) at all points in this plane except for a 
singular point at the origin where the system is metallic. 
Thus, this cyclic adiabatic change in H takes the system 
along an insulating path that encircles this singular point, 
so that a quantized particle transport AP — J(t) dt 
of a unit charge is obtained^ 

Away from the adiabatic regime, deviations from exact 
quantization are expected. This can be understood from 
the fact that under nonadiabatic conditions the state 
at time t depends on the history at times t' < t. In 
particular, the final state may be be different from the 
initial one even though H°(T) — H°(0), in which case 
P(T) — .P(O) 7^ 1. By contrast, an adiabatically-evolving 
system has no memory, being completely determined by 
the instantaneous H°(t) and dH°/dt. 

To illustrate this point, we increased a from to 2tt 
during a time interval t 6 [0, T] according to a(t) — 
2tt sin 2 (7rt/2T), and held it constant afterwards. The 
system was prepared at t = in its ground state, and 
the wave functions evolved in time according to Eq. (|(j5fl , 



FIG. 2: Time evolution of the polarization as a result of 
changing the sliding parameter a from to 2n over the time 
interval [0,80], using 200 k points. The solid line shows the 
actual dynamic polarization, while the dashed line shows the 
ground-state polarization of the instantaneous Hamiltonian. 
Inset: Detail of the remnant oscillations of the polarization 
after the Hamiltonian stops changing, at t = 80. For compar- 
ison, the results using 100 k points are also shown. 

using At = 0.005 (the same time step was used in all 
other simulations in this work). At each time step we 
computed the dynamic polarization using Eq. i|52[) . 

The resulting P(t) for T = 80 and 200 k points is 
shown in Fig. [3 where we also display the exact adi- 
abatic (T — > oo ) limit -P s tatic [<*(<)] obtained by diago- 
nalizing H®(a(t)) on the same mesh of k points^ (To 
check that our calculations are converged with respect to 
the number of k points, the inset of Fig. [21 compares the 
results for 100 and 200 k points.) The dynamic polariza- 
tion P(t) obtained by solving the TDSE follows closely, 
but not exactly, the adiabatic curve. In particular, at 
the time t — 80 when the Hamiltonian stops changing, 
the polarization differs slightly from unity (see inset in 
Fig. [2J) , indicating that the system is not in the ground 
state. The oscillations that follow arise from quantum in- 
terference (beats) between valence and conduction states, 
as a result of having excited electrons across the gap dur- 
ing [0, T]. Their period, of 5.5 time units, corresponds to 
the fundamental gap in Fig. ^ E gap = 1.137. This is 
consistent with the fc-space distribution of the (small) 
electron-hole pair amplitude present in the system after 
time T: for a sliding period of T = 80, the distribution 
is mostly concentrated around k = 0, and it is essen- 
tially the lowest conduction band that gets populated. 
As the adiabatic limit is approached by increasing T, the 
amplitude of the remnant oscillations of the polarization 
decreases. This is illustrated in the upper panel of Fig. [31 
where we compare T — 80 with T = 120. 

Besides the macroscopic polarization P(t), another 
quantity of interest is the electronic localization length 
£(t) that characterizes the root-mean-square quantum 
fluctuations of the macroscopic polarization^! It is given 



11 




0.4 0.6 
cc(t)/27t 



FIG. 3: Upper panel: Same as Fig. but now using two 
different time intervals, [0,80] and [0, 120], for changing the 
sliding parameter a. Lower panel: electron localization length 
versus the instantaneous value of a, during the time in- 
tervals [0, T] during which a is changing. 
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FIG. 4: Time evolution of the polarization as a result of in- 
creasing the electric field from to f max = 0.025 over two time 
intervals, [0, 40] and [0, 80], using 200 k points. The solid line 
shows the actual dynamic polarization, while the dashed line 
shows the static polarization for the instantaneous value of 
the field. Inset: Comparison of the dynamic polarization for 
100 and 200 k points. 



Gradual turn-on of an electric field 



by £ 2 = Clj/M, where fl\ is a gauge-invariant quantity 
which in one dimension is equal to the spread of the 
maximally- localized Wannier functions (|18[l . 3 We have 
computed fii using Eq. (34) from Ref . 0, and in the lower 
panel of Fig.|31we plot £(t) against a{t)/2i:. In the adia- 
batic limit the resulting curve consists of three identical 
oscillations, reflecting the existence of three equivalent 
atoms in the unit cell. As nonadiabaticity increases, £ 
tends to increase as well. Nevertheless, the electrons re- 
main localized, i.e., "insulating-like", in the sense dis- 
cussed in Sec. UnB^a 

The above results are representative of the regime 
where deviations from adiabaticity are small. If we in- 
crease the degree of nonadiabaticity by choosing a smaller 
T (e.g., T = 40), we begin to notice a linear increase of 
the polarization at later times. This new behavior can be 
traced to the excitation of electron and hole wavepackets 
centered at some ko and propagating at different group 
velocities. Let Ai?fc be the interband separation, and 
Av g be the difference of group velocities of the two low- 
est bands, at fco- Then, in addition to the quantum beats 
of period 2irh/ A£ , / Co caused by the interband dynamics, 
we observe a linear-in-t term in P(t) with slope propor- 
tional to Av fl , reflecting the change in dipole moment 
as the electron-hole pair separates. (More precisely, the 
preceeding statements apply only in the limit of a dense 
fc-point mesh; for any finite mesh spacing Afc, the linear 
behavior is replaced by an oscillatory one with an am- 
plitude scaling as 1/Afc and period 27r/(Aw s Afc). Thus, 
an especially fine fc-point mesh should be used if these 
effects are to be investigated.) 



In the previous example the electric field was held at 
zero, and the dynamics was produced by varying the 
parameter a in H . Let us now study the polariza- 
tion response of the system when an electric field £(t) 
is switched on linearly over a time interval [0, T] and is 
held fixed afterwards. We have set a = 0, so that the 
ground state is centrosymmetric, with zero spontaneous 
polarization. 

We begin by considering a situation where the final 
value of the field, £ ma xj is smaller than the /c-mesh- 
dependent critical field £ C (N) above which the energy 
functional l|45|l has no minima. This allows us to compare 
the dynamic polarization P(t) with the static polariza- 
tion P s tatic of the stationary state in the presence of 
the same field, which we find by minimizing the energy^ 
In Fig. 01 we display the results for £ max = 0.025 and two 
different switching times, T = 40 and T = 80. The simu- 
lation was done using 200 k points, to which corresponds 
a critical field £ C (N = 200) ~ 0.037. (The inset shows 
the agreement between the results obtained using 100 
and 200 k points.) Clearly, P(t) tracks quite closely the 
adiabatic curve -Pstatic [£(*)]) the more so as T increases. 
This illustrates the point, emphasized in Ref. |j, that the 
state obtained by minimizing a field-dependent energy 
functional should be thought of as the one which is gen- 
erated from the zero-field state by adiabatically turning 
on £. 

Let us now explore the regime above £ C (N), where 
energy- minimization schemes fail. For fmax > £ C (N) 
the exact adiabatic limit of the process of ramping up 
the field is unattainable. Nevertheless, if £ max is small 
compared to the field scale at which intrinsic breakdown 
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occurs (i.e., at which the Zener tunneling rate becomes 
on the of order interband frequencies, which is a bulk 
property^), a quasistationary state should be reachable 
by turning on the field at a rate that is slow compared 
to the usual electronic processes, but fast compared to 
the characteristic tunneling rate at the maximum field 
encountered. After the ramp-up is completed, but at 
times still short compared to the tunneling rate, this state 
should provide the appropriate extrapolation to fields 
above £ C (N) of the truly stationary state that exists be- 
low £ C (N). 

To illustrate this situation, we repeated the calcula- 
tion with 200 k points depicted in Fig. 0J but increasing 
the maximum field from 0.025 to 0.05, somewhat larger 
than £ c (200) - 0.037. The resulting curve for P(t) is 
very similar to that in Fig. without any sign of run- 
away behavior. As a more striking example, we show in 
Fig. [S] the outcome of calculations with the same final 
field of £ max = 0.05, but with even denser sets of 400 
and 800 k points. For the latter £ c ~ 0.01, considerably 
smaller than £ maX j and still there is no sign of instability. 
(Note also that the P(t) curve in Fig. [S] - whose verti- 
cal scale differs from that in Fig. 0] by the same factor 
of two that exists between the respective values of £ m ax 
- looks almost identical to that in Fig. 21) These results 
confirm that, as long as we are solving a time- dependent 
Schrodinger equation for a given history of switching on 
the field, there is no such thing as a Afc-dependent criti- 
cal field; the thermodynamic limit of an infinitely dense 
fc-point mesh is perfectly well defined. The only break- 
down behavior that may be observed in short time scales 
is the physical one that occurs when the applied field 
is large enough that the Zener tunneling rate becomes 
significant 42^i The concept of a Afc-dependent critical 
field applies only to the attemp to obtain solutions in the 
presence of a static electric field from an energy varia- 
tional principle. By going back to the original dynamical 
problem of slowly ramping up the field, we circumvent the 
difficulties that ultimately resulted from trying to treat 
as a (stable) stationary state what is really a long-lived 
resonance. 



D. Dielectric function in a static field 

There is great interest in modulating the optical prop- 
erties of crystals and superlattices by applying static 
electric fields. An example of such an electro-optical ef- 
fect is the modification of the dielectric function. This 
is known as the Franz-Keldysh effect, or electroabsorp- 
tion. Although it has been extensively studied in bulk 
semiconductors^ quantum wells^ and superlattices^ 
we are not aware of any first-principles investigations. 
The present method may provide a route to such calcu- 
lations. 

We compute the dielectric function in the presence of a 
static field £o as follows. The system is prepared at t = 
in the stationary state polarized by a field £q + A£ , with 
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FIG. 5: Same as Fig. but now using £ max = 0.05 and 800 
k points. Since £ max is larger than the critical field for this 
number of k points [S c ~ 0.01), no adiabatic curve Pstatic [S(t)] 
is shown. Inset: Comparison of the dynamic polarization for 
400 and 800 k points. 
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FIG. 6: Susceptibility x [£ol H in 

the presence of a static field 
£o, for a = and 100 k points, using a level-broadening S = 
0.04. Dotted lines: Kubo formula result for So — 0; solid lines: 
results using our method, for both £ o = and £o = 0.05. The 
latter displays the Franz-Keldysh effect. The inset compares 
the Franz-Keldysh oscillations for two different bias fields, 
So = 0.05 and £ = 0.03. 
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\A£\ << \£o\. By using a field of magnitude below the 
critical field, we are able to find that state by minimizing 
the energy. For t > we let the system evolve in time 
in the presence of the target field So- Let -P s tatic[£o] be 
static polarization of the system under the field £q. The 
polarization response to the step-function discontinuity 
in £{t) = £ Q + A£9{-t) is AP(t) = P(t) - P 8tat ic[£o]- 
To obtain the frequency-dependent response we need the 
Fourier transform of AP(t) for t > only: 



r+oc 

AP{u) = / AP(t) e^-^ 1 dt, 
Jo 



(75) 



where a damping factor S has been introduced as an ap- 
proximate way to account for level broadening. 45 To lin- 
ear order in A£ the susceptibility is 



Rex [£ol H 



dP Ki 



d£ 



£=£o 



A£ 



ImAP(w), (76) 



Imx [£ol M 



A£ 



ReAP(w) 



(77) 



With this real-time approach the need to perform a 
summation over conduction-band states is circumvented. 
Previous real-time, scalar potential approaches" were 
restricted to finite systems, since it was unclear how 
to evaluate the dynamic macroscopic polarization of an 
extended system. A real-time, vector-potential scheme 
valid for bulk systems was proposed in Ref. ItH 

We validate our method by comparing in Fig. 
the ground-state susceptibility with the analytic Kubo- 
formula (sum-over-states) result, using in both cases the 
same broadening 5 and fc-point mesh. Also shown in 
Fig. Elis the susceptibility in the presence of a £q = 0.05 
bias field, displaying the Franz-Keldysh effect: an absorp- 
tion tail below the gap caused by photon-assisted tunnel- 
ing, and oscillations above the gap42. The Franz-Keldysh 
oscillations become more widely spaced with increasing 
£q. This is illustrated in the inset of Fig. where we 
compare them for £q = 0.05 and £q = 0.03. 



VII. SUMMARY 

The work of King-Smith and Vanderbilt demonstrated 
that the bulk electronic polarization, defined in terms 
of the current flowing during the adiabatic evolution of 
an insulating system in a vanishing macroscopic elec- 
tric field, could be related to a Berry's phase defined 
over the manifold of occupied Bloch statesii We have 
generalized this result by considering the time evolution 
of an initially insulating electron system under the very 
general Hamiltonian (J2J), where the lattice-periodic part 
H°(t) and the homogeneous electric field S(t) may have 
an arbitrarily strong and rapid variation in time. In 
the absence of scattering, we have proved that the in- 
tegrated current AP = J J(t) dt is still given by the 



King- Smith- Vanderbilt formula, but written in terms 
of the instantaneous Bloch-like solutions of the time- 
dependent Schrodinger equation. The coherent dynamic 
polarization P(t) was interpreted as a nonadiabatic ge- 
ometric phased These generalizations of the theory al- 
lowed us to justify recent developments in which the en- 
ergy functional E of Nunes and Gonzei has been used 
as the basis for direct DFT calculations of insulators in 
a static homogeneous electric fieldAS The limitation of 
those methods to fields of magnitude smaller than a Ak- 
dependent critical field that vanishes in the thermody- 
namic limit has been removed: we have shown numeri- 
cally that quasistationary states in finite fields exist for 
arbitrarily dense fc-point meshes, and can be obtained 
by solving the time-dependent Shrodinger equation for a 
slowly-increasing field. The present method also provides 
a convenient framework for the computation of coherent 
time-dependent excitations in insulators. As an example, 
the dielectric function was calculated for a tight-binding 
model by considering the response to a step-function dis- 
continuity in £(t), illustrating effects such as photon- 
assisted tunneling and Franz-Keldysh oscillations. A full 
ah initio implementation within the framework of time- 
dependent density-functional theory should be possible. 
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APPENDIX A: CRYSTAL-MOMENTUM 
REPRESENTATION 

The introduction of linear scalar potentials in crys- 
tals is usually discussed in the language of the crystal- 
momentum representation (CMR). 14 Instead, we have 
used the Berry-phase theory of polarization, and the pur- 
pose of this Appendix is to show how to switch from one 
to the other. The CMR uses as a basis the eigenstates 
IV'km) of H° with eigenvalues E^ m . In accordance with 
Eq. (fP7|) we assume that IVWiC 1 *)! 2 integrates to one over 
the unit cell volume v. That implies 48-49 



(ipkm\ipk'i) = J ^Zmi^uairfdr = n B S(k - k')S m i. 

(Al) 

The CMR expansion of the identity operator is 

i = fi- 1 ^ / dk|Vw><Vw|, (A2) 

m—l 

so that a general one-electron state \(f>) is expanded as 

oo „ 

|0) = i\4>) = Y, j dk 7k' m |^k' m ), (A3) 



14 



where /k' m = ^B 1 (V , k'm|</')- For the occupied Bloch-like 
states \4>kn) m a WR manifold, the CMR wave function 
/k'm(k, n) takes the form 

/k' TO (k, n) = c k , inm S(k' - k) (A4) 

with Y^m=l l c k,nm| 2 = 1j which leads to Eq. ©. 

1. Current and the CMR velocity operator 

The velocity operator l|31|l is diagonal in k and is con- 
veniently split into a sum of two operators, one diagonal 
and the other off-diagonal in the band indexsi^ 



v = v d + v od . 



(A5) 



The matrix elements of v d are 

(Vw|v d |Vk<z) = n B 6(k - k')<5 mi v d m , (A6) 

where 



1 



<9k^kr, 



(A7) 



The matrix elements of v od are 

(^ m \v° d \M = n B 5(k - k')< mi , (as) 

where 



V k d mJ = ^ X k,m2 [Ekm - E k i] 



(A9) 



and we have defined the hcrmitian matrix 

x u, m i = i (ukm\d ka u k i) , (A10) 

which is analogous to Eq. 1)26(1 for the |i>kri)- 

The current, Eq. (|28|l . is split into intraband and in- 
tcrband parts, 

inter 

(*)• (AH) 

Writing the density matrix as 

(ipkm\n\ipk'i) = ^B<5(k - k')n k ,mz, (A12) 

where 

M 



^k,mZ ^ ^-k.nm [^k,n/] ; 



(A13) 



n=l 



we find 



Jintra = Tr c (nv d ) = ^ / dkn w . mm v^ 

v ^ m=l ^ 



(A14) 



and 

Jintcr 



^Tr c (nv° d ) = --^ £ / «ikn k , mi v k 



d 
k,lr 



In the above we used the CMR form of Eq. i(2"§|l , 

Tr c (0) = 03 1 ^ f dk^(^ m \0\i;u m ), (A16) 

m=l ^ 

where iV should be taken to signify UbS(0)^ 

Plugging into l(27|l yields, after some manipula- 
tions, Eqs. (|A11(I . (|A14|) . and (|A15|) . confirming that the 
Berry-phase polarization correctly accounts for both in- 
traband and interband contributions. It is instructive 
to consider some particular cases. The adiabatic current 
J = (dP/dX)\ discussed in Refs. flll3^ is purely interband. 
If the perturbation is a sinusoidal electric field, the linear 
response is again a purely interband current, while the 
nonlinear response has also an intraband component^S^i 



2. Polarization and the CMR position operator 

Along the same lines, one can show that the Berry- 
phase expression for P is consistent with the CMR posi- 
tion operator, which takes the formM 

(V>km|*#k'/) = -«^B<9k'<5(k' -k)S mn +n B 8(k' -k)X k , m ;. 

(A17) 

Combined with Eqs. (|A12|) and (|A16|) this yields 



00 „ 

P = - -Tr, (;/r) = ^3 J rfkn k , m; X k 



(A15) 



v ■ ■ (2tt) 3 ^ 

which is the same results one gets from inserting the 
CMR expansion 10 into the nonadiabatic Berry-phase 
formula (|21ll . The linear character of f is reflected in the 
above equation being defined only up to a quantum of 
polarization. 



3. CMR dynamical equations 

In the case where H° (and hence the CMR basis) is 
constant in time, plugging © into the TDSE yields 
the CMR form of the Schrodinger equationj22i£2 

iTic^m = (E km + ie£ ■ V k )ck m , (A19) 

where we ticive simplified Ck,nm 

to Ckm and defined 

^kCkm = <9 k C k m - i ^2 X k;Ck/, (A20) 

1=1 

which is reminiscent of the covariant derivative, Eq. 142|) 
(but note the difference in the sign of the last term). It 
is customary to write Eq. (|A19|) as 

00 

ihc km = (E^+ieS-dkjc^jn+eS ■ ^ C k/X k , m z, (A21) 

1^771 
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where 



(A22) 



is a shifted energy eigenvalue. E^ is identical to Eq. lpT7|) 
except that i>km) has been replaced by the zero-field 
eigenstate Ukm). Upon averaging over k the last term on 
the right-hand-side becomes the first-order shift in total 
energy, —vPo ■ £, where Po is the spontaneous Berry- 
phase polarization. 

In general the above TDSE has no stationary solutions. 
Approximate solutions - the Wannier-Stark states - re- 
sult from restricting the wavepacket dynamics to a single 
band (the semiclassical approximation) . That is achieved 
by dropping the sum on the right-hand-side of Eq. I|A21|) . 
which is responsible for interband tunneling4ii£2 

Finally, combining (|A13|I and l|A19|) produces the dy- 
namical equation for the CMR density matrix: 

^ %,nm \Ekn ^km)^k,nm ~t~ 2e£ • t^k^k, nm 

oo 

— e£ ■ («k,niXk,im — Xk,„znk,im)- 

1=1 

(A23) 

A closely related form has been used to study the non- 
linear optical susceptibilities of semiconductors , 5 1 1 53 



APPENDIX B: COVARIANT DERIVATIVE AND 
RELATED OPERATORS 



In Sec. lIV Al we introduced a modified TDSE that con- 
tains the multiband covariant derivative dk, Eq. (|42|l . 
that was instrumental for making contact with the 
discrete- k dynamical equations of Sec. 11V Bl Here we 
summarize the properties of the covariant derivative and 
other closely related operators. 

The covariant derivative <9k|«kn) of an occupied state 
transforms in the same way as that state under a gauge 
transformation, Eq. (JTDJ: 



(Bl) 



Moreover, it is orthogonal to the occupied subspace at k, 



(Vkm\dkV kn ) = 0. 



(B2) 



Recalling that parallel transport is characterized by 
(vkn\dkVkn) = 0, for m = n this relation shows that 
dk acting in an arbitrary gauge gives the same result as 
(?k acting in the parallel-transport gauge that shares the 
same states at k. In the discretized form (|(j(jp) the prop- 
erty I|B1|I is a consequence of Eq. lf531) , and the property 
(IB2II is a consequence of Eq. (|54f) . Like i&k, i&k is her- 
mitian. By this we mean that its matrix representation 



in an orthonormal basis (e.g., the |«kn), n = 1, ...,M 
complemented by a set of unoccupied states |ckj)) is her- 
mitian. This is closely related to the hermiticity of the 
matrix A£ defined in Eq. I|26|) . Finally, note that 



i<9k|«kn) = iQkdk\vkn) = iQAPklvkn) , 



(B3) 



i.e., the action of id^ on an occupied state is identical 
to that of iQk<9k and iQ^d^P^. They differ in how they 
act on the unoccupied states. Unlike idk, the other two 
are not hermitian: for instance, («Qk<9kAc)^ \vkn) = 0. It 
follows from these considerations that Eq. (|41|) can be 
recast as 

ih\v^ n ) = [Hi + e£ ■ {iQ k d k P k + h.c.)] |« kn ). (B4) 

This is the form of the TDSE to which Eqs. 
reduce in the continuum- fc limit, since 



w k ~ ie£ ■ Qk<9 k -Pk 
(compare with Eq. H59|) ~). 



(B5) 



APPENDIX C: GRADIENT OF THE ENERGY 
FUNCTIONAL 

The purpose of this Appendix is to obtain expressions 
for the derivatives of the two terms in the energy func- 
tional of Eq. (1^31) with respect to the occupied Bloch-likc 
states in the discrete-fc case. The results have been used 
in Sees. IIVBI and IV Al for the discussion of the time- 
dependent evolution equations and the stationary solu- 
tions, respectively. 



1. Band-structure contribution 

To find the gradient 5E/{5v-k n \ of the energy functional 
(|45|l . let us isolate the terms that depend on (wknl- Using 
(|43|l the zero-field part l|46l) can be expressed as E = 
(1/AQ E k so that 



5E° _ 1 <Str[P k #£] 



(Ql) 

{6v kn \ N (5v kn \ 

In order to allow for arbitrary variations of (wknli even 
those for which the (i>kn | do not remain orthonormal, we 
write 

M 

P*= E iS^) mn \vk m ){v^ n \, (C2) 

m,n— 1 

where 5k, mn = (vkm\vkn)- Dropping the subscript k, 
Str[PH°] = tr[(5P)H°] 

= E {S~ l ) mn [(v n \H°\Sv m ) + {Sv n \H°\v m ) 



-J2(v n \H°\v m )5(S-% 



(C3) 
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Using 8(S x ) = —S 2 SS and 8S mn — {v m \8v n ) 
(8v m \v n ), and evaluating at S = 1, we arrive at 



<5tr[P k g£] 



(C4) 



Thus the consequence of expressing Pk as (|C2() instead 
of (|43JI is to render the gradient orthogonal to the oc- 
cupied manifold at k. (When we derived the dynamical 
equation (|58|l using (|47|l , the gradient of E° was not or- 
thogonalized, which is why the dynamics did not follow 
parallel transport (see Sec. IIV (ty ). 



The corresponding derivative of the last term of l|C9p van- 
ishes since 5(k',k) = (^k'nl^kn) does not contain (i>kn| 
as a bra. We thus arrive at 



which combined with (|C7ll gives 



(C12) 



2. Polarization contribution 

To find the gradient of the field-coupling term — vS ■ P 
we need STi j (Sv^n | ■ Let us start by recasting Eq. I|52|) as 



8T % 



k« k« 



1 = «T X! X! ' ' "/ 

4 (=1 j=0 

where we have defined the phase 

^(k,k') = -ImlndetS*(k,k') 
Using 0(k',k) = — </>(k, k'), this becomes 



Aki), 



(C5) 



(C6) 



(C7) 



where only the terms depending on (i>kn| were written 
explicitly. Hence 



sr,. 



AT-i- 51 



k, kicr). 



(C8) 



The phase </>(k, k') can be expressed as 
<Kk,k') = -Imtrln5(k,k') 

= -trm5(k,k') - -trln5(k',k). (C9) 

For an arbitrary non-singular matrix A we have 

<5tr InA = tr \n(A + 8 A) - tr ln(A) 

= tr ln[(A + 8 A) A- 1 ] = tr ln[l + (8 A) A' 1 } 
= ix[A- 1 8A] + 0(5A 2 ), (CIO) 



so that 



gtajngfo k' 



= tr 

M 



^(k,kO^ (k ' k/) 



(^kn| 

= 5 mn( k ' k ') l^k'm) 



m— 1 

ISifn). (CH) 



(^ k „| 2^^ 



c|lW,n)- 



(C13) 



This is automatically orthogonal to the occupied man- 
ifold at k. (See Refs. [7HT0I for alternative derivations.) 
Collecting terms and using Eq. I|57[) . we obtain Eq. I|69|) 
for the gradient of the full energy functional E. 



APPENDIX D: DISCRETIZED FORMULA FOR 
THE CURRENT 



Just as the macroscopic polarization P is evaluated in 
practice via a finite-difference formula on a mesh of k 
points, the same can be done for the macroscopic cur- 
rent J = dP/dt. The invariance of Eq. (|27|l under the 
replacement <9k — > allows us to then use the discretiza- 
tion rule l|t)U|) . leading to 



J = ^EE H ^K»|-ffkK«,,n>ai + C.C., 
n,k i—1 <y— ±1 * 

( D1 ) 

which is significantly easier and cheaper to compute 
than the spatial average of the microscopic current. Wc 
have checked numerically on our one-dimensional tight- 
binding model that Eq. (|Dlf) yields, for small At, the 
same result as [P(t + At) — P(t)]/ At computed with the 
discrctizcd Berry-phase formula. 

The same strategy as outlined above can be used 
to derive a discretized formula for the Berry curvature 
(|25|l summed over bands, which is also invariant under 
3k — > 9k. This may be useful in other contexts, such as 
semiclassical wavepacket dynamics in crystals^* 
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